Explore genetic evolution only:

Constant Environment

To explore this we used the following parameters space :

constant_env=getAlllSummaries("unifPop50k")

I was not keeping track of the outpute rate so the convertion in real generation has to be done by hand:

constant_env$outputrate=10
constant_env=updateScale(constant_env)

Check the trajectories of the simulation for the population size for different mutation rate \(\mu\)(from left to right: \(\mu = 0, \mu = 0.001,\mu=0.01\)), different selective pressure \(\sigma_s\) (from top to bottomright: \(\sigma_s = 0.1, \sigma_smu = 0.2,\sigma_s=0.4,\sigma_s=1000\))and different carrying capacity \(K\)

plotAllTrajVar(constant_env,m=0.05,E=0,obs="N")

Check the variance at the

plotAllVariableSummaries(constant_env,E=0,estimate=eq2833b) #we use only simulations without noise

Constant but noisy environment

To explore this we used the following parameters space :

v=getAlllSummaries("growingDelta",exclude="outputrate")
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
noisy_env = rbind(v,constant_env[,colnames(v)])
noisy_env = noisy_env[ noisy_env$mu >0,]

To compare it we merge the two tables of results, and remove simualtion with no mutations

Check the trajectory off the effective population size:

plotAllTrajVar(noisy_env,m=0.05,E=0,obs="var_x")

Check the trajectory off the mean variance size:

plotAllTrajVar(noisy_env,m=0.05,E=0,obs="var_x",ylim=c(0,.005))

Check and compare the final variance :

plotAllVariableSummaries(noisy_env,E=0,estimate=eq2833b) 

Some of the data seems to be missing because when \(\delta>\sigma\) we have extintions, has shown in the graph:

sigmas=unique(noisy_env$sigma)
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
tmpc=noisy_env[ noisy_env$sigma==0.4 & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == 0))
for(s in sigmas){
    tmpc=noisy_env[ noisy_env$sigma==s & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
    Ks=sort(unique(tmpc$K))
    for(k in 1:length(Ks)){
        t=tmpc[tmpc$K==Ks[k],]
        r=tapply(t$extinction,t$delta,mean)
        lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
    }
}
        legend("right",
               legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
               col=c(rep(1,length(Ks)),cols),
               lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
               pch=c(seq_along(Ks),rep(NA,length(sigmas)))
               )

#dev.off()

Thus we introduce autocorrelation in the environment. To introduce autocorrelation we increase slightly \(omega\)

par(mfrow=c(1,3))
plot(1:500,environment(500,omega=1,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=2,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=3,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))

Constant, noisy, with autocorrelation

For this we re-run previous experiment but with \(\omega \in {1,2}\)

omega1=getAlllSummaries(c("omega1growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega1=updateScale(omega1)
omega1$omega=1

omega2=getAlllSummaries(c("omega2growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega2=updateScale(omega2)
omega2$omega=2

noisy_env_omega = rbind(v,constant_env[,colnames(v)],omega1[,colnames(v)],omega2[,colnames(v)])

Let’s look at the extinction again for the different omegas

par(mfrow=c(1,3))
omegas=unique(noisy_env_omega$omega)
sigmas=sort(unique(noisy_env_omega$sigma))
deltas=sort(unique(noisy_env_omega$delta))
Ks=sort(unique(noisy_env_omega$K))
for(o in omegas){
    cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
    names(cols)=sigmas
    tmpc=noisy_env_omega[ noisy_env_omega$sigma==0.4 & noisy_env_omega$mu == 0.01 & noisy_env_omega$omega == o & noisy_env_omega$m == 0.2 ,]
    xrange=range(noisy_env_omega$delta[noisy_env_omega$omega>0])
    plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == .(o)),xlim=xrange)
    for(s in sigmas){
        tmpc=noisy_env_omega[ noisy_env_omega$sigma==s & noisy_env_omega$mu == 0.01 & noisy_env_omega$m == 0.2 & noisy_env_omega$omega == o,]
        for(k in 1:length(Ks)){
            t=tmpc[tmpc$K==Ks[k],]
            r=tapply(t$extinction,t$delta,mean)
            lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
        }
    }
    legend("bottomleft",
           legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
           col=c(rep(1,length(Ks)),cols),
           lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
           pch=c(seq_along(Ks),rep(NA,length(sigmas)))
           )
}

We can thus reproduce he trajectories and variance equilibrium using simulation with \(\omega=1\) et \(\omega=2\)

Check the trajectory off the effective population size:

#pdf("trajectory_N_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="N")

#dev.off()
#pdf("trajectory_N_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="N")

#dev.off()

Check the trajectory off the mean variance size:

#pdf("trajectory_var_x_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="var_x",ylim=c(0,.002))

#dev.off()
#pdf("trajectory_var_x_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="var_x",ylim=c(0,.002))

#dev.off()

And check the variance at the end of the equilibrium:

constant_env$omega=0
omega1_wdz=rbind(omega1,constant_env[constant_env$mu>0,])
omega2_wdz=rbind(omega2,constant_env[constant_env$mu>0,])
omega1_wdz=omega1_wdz[omega1_wdz$sigma<10,]
omega2_wdz=omega2_wdz[omega2_wdz$sigma<10,]
plotAllVariableSummaries(omega1_wdz,E=0,estimate=eq2833b) #we use only simulations without noise

plotAllVariableSummaries(omega2_wdz,E=0,estimate=eq2833b) #we use only simulations without noise

Moving optimum \(\theta_t=vt\)

par(mfrow=c(1,3))
vts=c(0.001,0.002,0.003)
cols=colorRampPalette(c("darkgreen","green"))(length(vts))
names(cols)=vts
omegas=1:3

for(o  in omegas){
plot(1,1,xlim=c(1,500),ylim=c(-.1,1.5),type="n",xlab="t",ylab=expression(theta[t]),main=bquote(omega == .(o)))
     for ( vt in vts)
         lines(1:500,environment(500,omega=o,delta=.1,vt=vt),col=cols[as.character(vt)])
legend("topleft",legend=paste0("v=",vts),col=cols,lty=1)
}

moving_theta=getAlllSummaries("growingThetat")
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
moving_theta=updateScale(moving_theta)
par(mfrow=c(1,4))
rates=unique(moving_theta$vt)
sigmas=sort(unique(moving_theta$sigma))
deltas=sort(unique(moving_theta$delta))
Ks=sort(unique(moving_theta$K))
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
    names(cols)=sigmas
for(d in deltas){
    tmpc=moving_theta[ moving_theta$sigma==0.4 & moving_theta$mu == 0.01 & moving_theta$delta == d & moving_theta$m == 0.2 ,]
    #xrange=range(moving_theta$delta)
    plot(tmpc$vt,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == .(o)))
    for(s in sigmas){
        tmpc=moving_theta[ moving_theta$sigma==s & moving_theta$mu == 0.001 & moving_theta$m == 0.2 & moving_theta$deltas == d,]
        for(k in 1:length(Ks)){
            t=tmpc[tmpc$K==Ks[k],]
            r=tapply(t$extinction,t$vt,mean)
            lines(sort(unique(t$vt)),r,type="b",pch=k,col=cols[as.character(s)])
        }
    }
    legend("bottomleft",
           legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
           col=c(rep(1,length(Ks)),cols),
           lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
           pch=c(seq_along(Ks),rep(NA,length(sigmas)))
           )
}